Phase diagram of the toric code model in a parallel magnetic field 
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Ground-state phase diagram of the toric code model in a parallel magnetic field has three dis- 
tinct phases: topological, charge-condensed, and vortex-condensed states. To study it we consider 
an implicit local order parameter characterizing the transition between the topological and charge- 
condensed phases, and sample it using continuous-time Monte Carlo simulations. The corresponding 
second-order transition line is obtained by finite-size scaling analysis of this order parameter. Sym- 
metry breaking between charges and vortices along the first-order transition line is also observed, 
and our numerical result shows that the end point of the first-order transition line is located at 
hi c) =h{ c) =0.418(2). 
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I. INTRODUCTION 

Topological phases are attracting a lot of attention 
for a variety of reasons from precise determination of 
physical constants to unusual properties of quasiparti- 
cles and promise for quantum information and computa- 
tion applications. The search for such phases has been 
mainly focused on electron systems, such as quantum 
Hall states and topological insulators 1 ' 2 . However, topo- 
logical phases can also exist in magnetic systems and the 
toric code model (TCM) ! based on spin- 1/2 degrees of 
freedom provides a key example. TCM can be viewed 
as the anisotropic limit of Kitaev's honeycomb model', 
which has recently been suggested to emerge in iridium 
oxides 5- ' . 

The ground state of TCM in zero external field is 
topologically ordered and characterized by two types of 
bosonic excitations (charges and vortices) with mutual 
scmionic statistics. To be useful for quantum comput- 
ing, the topological character of TCM should be robust 
against local magnetic field perturbations. Thus one of 
the most important questions to ask is what are the criti- 
cal values of magnetic fields which cause condensation of 
charges or vortices, i.e. break the topological order. The 
purpose of this paper is to perform an accurate study 
of phase transition lines in TCM subject to a parallel 
magnetic field. 

As TCM involves four-spin interactions, the 
continuous-time Monte Carlo approach of Ref. 8 
has to be modified appropriately and supplemented with 
additional updates. The topological phase is charac- 
terized by non-local quantum correlations. However, 
an implicit local order parameter can still be defined 
to characterize the transition between the topological 
phase and charge-condensed phase. The corresponding 
second-order transition line is then obtained from the 
standard finite-size scaling analysis of this order param- 
eter. Along the self-duality line, there is a symmetry 
between charges and vortices which is, however, broken 
in a certain range of parameters. It is observed in 



the Monte Carlo simulation as the first-order phase 
transition between the charge- and vortex-condensed 
phases. 

Similar studies of TCM were performed in the past us- 
ing other methods' -12 . In Ref. 9 the authors used an 
approximate mapping of TCM onto the anisotropic Z2 
gauge Higgs model, and computed the phase diagram 
of the corresponding 3D classical model in the isotropic 
limit, extending previous studies of the Z2 gauge Higgs 
model, see Ref. 13, to larger system sizes. The behav- 
ior of TCM in a parallel field and/or a transverse field 
was also studied using several perturbative (high order) 
treatments 1 " -12 . The Monte Carlo algorithm developed 
in this paper suffers from the sign-problem in the pres- 
ence of the transverse field component. Therefore we 
restrict our studies to the parallel field in which case 
the parameter space of the model considered here is the 
same as in Ref. 10. Not surprisingly, we find the same 
topology of the phase diagram, but our numerical result 
shows that the end point of the first-order transition line 
is located at h { x c) = h[ c) = 0.418(2), which differs from 

hi c} = h { z c} = 0.48(2) calculated in Ref. 10. To elabo- 
rately compare our Monte Carlo data with the series ex- 
pansion results, we also compute two magnetization-like 
quantities using the bare series expansion given in Ref. 10 
and 12. We argue that the discrepancy in the end point 
is mainly because the series expansion about high-field 
limit at order 5 is questionable in the field strength of 
interest. 



II. MODEL 

The toric code model is defined on a square lattice with 
spins located on lattice bonds. The Hamiltonian (in zero 
field) is based on four-spin interactions: 
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FIG. 1: (Color online) An illustration of the toric code model 
geometry. A s = Ylj £s &j and B p — Ylj ep a j are products of 
spin operators on bonds incident to site s and surrounding 
plaquette p, respectively. 

Pauli matrices). Here subscripts s and p refer to sites 
and plaqucttes of the square lattice, see Fig. 1, respec- 
tively, and address all spins surrounding the correspond- 
ing site or plaquette. Since all A s and B p commute with 
each other, the ground state manifold is known exactly. 
It corresponds to eigenstates of all A s and B p with max- 
imal eigenvalues +1. The ground-state degeneracy de- 
pends on boundary conditions, and on a torus there are 
four degenerate ground states. Elementary excitations, 
a charge on site s or a vortex on plaquette p, correspond 
to states with eigenvalues A s = — 1 or B p = — 1. When 
a charge (vortex) is circulated around another vortex 
(charge) the wave function acquires a phase shift of ir. 
Thus charges and vortices have mutual semionic statis- 
tics. On the other hand, composite charge-vortex parti- 
cles have fermionic statistics. 

In the presence of a parallel magnetic field, which is 
perpendicular to the y— direction, the Hamiltonian takes 
the form: 

H = H TCM -KY, °l -KYs a l > ( 2 ) 

b b 

where b refers to lattice bonds. By orienting individ- 
ual spins the h z and h x terms reduce energy gaps for 
creating charge and vortex excitations, respectively, and 
ultimately cause their condensation. 



In the basis of of eigenstates, U and K are conveniently 
called potential and kinetic energy terms, respectively. 
Similarly to the standard Fcynman's lattice path integral 
in the interaction representation we expand the partition 
function Z in powers of K as follows: 

Z = Tr[exp(-/3iT)] 

oo „ 

= So<ti<T3<—<t„ ^ Tl " ' S{a ,ai— a n _i,a„=a } 
n— 

(-l)™^^ • • • K a2ai K aiao exp{- J* U{r)dT}. (4) 

Here ctj with j = 1,2 ■■■n enumerate of basis states 
characterized by a set of classical variables Zb = ±1 on 
all bonds; Tj are the imaginary time moments when the 
state Oj_i is changed to aj due to action of the kinetic 
energy term, K a . a . 1 = (aj \ K \ a.j-\). By definition, 
the four-spin operator A s simultaneously flips four spins 
on the bonds connected to the site s, while operator 
flips only the spin on bond b. Finally, the dependence 
of potential energy on time in the exponential of (4) is 
defined as U(t) = U ai for r € (r^Tj+i). Note that the 
weight for each configuration is positive here, while the 
inclusion of the transverse field component will result in 
the sign problem. 

The continuous-time Monte Carlo has been mainly 
used for quantum models with two-body interactions, 
such as the Bose-Hubbard model, where the worm 
algorithm 8 provides an effective approach to large-scale 
and high-efficiency simulations. Here we apply the 
method to quantum models with four-spin non-diagonal 
interactions. An obvious set of elementary Monte Carlo 
updates, which (1) create/annihilate a pair of one- 
spin flips, (2) change times of one-spin hips, (3) cre- 
ate/annihilate a pair of four-spin flips, and (4) change 
times of four-spin flips, is not ergodic for Hamiltonian (2) 
because it does not admit all allowed configurations with 
alternating single-spin and four-spin flips. This problem 
is taken care of by an update which involves both four- 
spin and one-spin flips. In the following, we provide a 
detailed description of all updates. 




III. CONTINUOUS-TIME MONTE CARLO 
SIMULATION 

To achieve complete quantum mechanical solution of 
the model we use the continuous-time Monte Carlo ap- 
proach introduced in Ref. 8. Within this approach, the 
Hamiltonian(2) is split into two non-commuting terms 
H = U + K with 

U = -J P YB P -h z J2<rb, 

P b , . 

(3) 

K= -J s Y A s -h x Y a b- 

s b 



FIG. 2: (Color online) Create a pair of one-spin flips. 

Update 1 Create/ Annihilate a pair of one-spin flips. 
Step 1.1 Select an bond b at random. Let N(b) denote 
the number of one-spin flips on bond b. 
Step 1.2 Keep of (r = 0) unchanged, and propose to cre- 
ate or annihilate a pair of one-spin flips on bond b with 
equal probability. 

(T) For creation, draw a random integer j = [(N(b) + 1) x 
rndm] (j = 0, 1 • • • N(b)) to select the time interval to be 
updated (rj b) , rj^), where with .] = !■■■ N(b) is the 
time position of the jth one-spin flip on bond b and = 
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0, 7jyLv +1 = P- Draw times for the new one-spin flips us- 
ing simple uniform probability density within the selected 



interval, W(r a ,Tp > r Q ) = const = 2/(rj+ ) 1 -rf'f. The 
acceptance ratio for creation is given by 
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A[/(r)dr}}, 



where AC/ (r) is the change of the potential energy (the 
same notation will be used for other updates without ex- 
plicitly mentioning that the integral has to be taken only 
for the updated time interval). The proposed configura- 
tion change is illustrated in Fig. 2. 



© In creation, draw a random integer j — [(N(s) + 
1) x rndm] (j = 0, l---N(s)) to identify an interval 
(rj s \ rj^), where Tj is the time position of the jth 

four-spin flip at site s with Tq = and Tj^) s \ + i = P- 
Propose time positions for the new four-spin flips from 
a simple uniform probability density within the selected 
interval, W(r a ,rp > r a ) = const = 2/(rj^ ) 1 -rj s) ) 2 . The 
acceptance is then 

f T W _ T ( S ))2 r 
fl = min{l, V 3+1 3 J 2 exp{- / AJ7(r)dr}}, 

Fig. 4 provides an illustration of the configuration change 
implied by the proposed update. 
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FIG. 3: (Color online) Annihilate a pair of one-spin flips. 

(2) For annihilation, exit the update if N(b) < 2; other- 
wise draw a random integer j = [(N(b) — I) x rndm] + I 
(j = 1 , 2 ■ • ■ N(b) — I) and propose to annihilate jth and 
(j + I)th one-spin flips. The acceptance ratio is 



R = 



(r 



(6) 
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^exp{-y AC/(r)dr}} 



The update is illustrated in Fig. 3. 
Update 2 Move a one-spin flip. 

Step 2.1 Randomly select a bond b. Exit this update if 
N(b) = 0. 

Step 2.2 Randomly select a one-spin flip, j = [N(b) x 
rndm] + 1 (j = 1, 2 • • ■ N(b)). Let its time position be t 
and times of previous and next flips on bond b (regardless 
of their types) be t\ and i 2 respectively, with noting that 
the worldline is 0— periodic. 
Step 2.3 

® If t\ < t < t-x, propose to move the flip to t' = (£2 — 
ti) x rndm+ti; 

(2) Otherwise, the new time position of the flip can be 
either in (0,t 2 ) or (ti,(3), i.e. t' = mod{(/3 + t 2 - h) x 
rndm + ti, /?}. 

Step 2.4 The acceptance ratio is simply given by the 
potential energy change 

i? = min{l,exp{- J MJ(t)<It}}. 



Update 3 Create/ Annihilate a pair of four-spin flips. 

This update is similar to Update 1. 

Step 3.1 Randomly select a site s. Let the number of 

the existing four-spin flips at s is N(s). 

Step 3.2 With a^(r = 0) on bonds b\, . . . , 64 attached 

to s (sec Fig. 4) kept unchanged, propose to create or 

annihilate a pair of four-spin flips at site s with equal 

probability. 
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FIG. 4: (Color online) An illustration of creating a pair of 
four-spin flips. 

(2) In annihilation, exit the update if N(s) < 2; otherwise 
draw a random integer j = [(N(s) — 1) x rndm] + I ( j = 
1, 2 ■ • • N(s) — 1), and propose to annihilate the jth and 
(j + l)th four-spin flips. The acceptance ratio is 

R = min{f , 2 j,exp{- j AU(r)dr}}, 

\ T j+2 T j-l) 8 

and the configuration changes are illustrated in Fig. 5. 



i+i 



FIG. 5: (Color online) An illustration of annihilating a pair 
of four-spin flips. 



Update 4 Move a four-spin flip. This update is similar 
to Update 2. 

Step 4.1 Select a site s at random. Exit this update if 
N(s) = 0. 

Step 4.2 Randomly select a four-spin flip on site s, j = 
[N(s) xrndm] + I (j = 1, 2 • • • N(s)). Let its time position 
be t while t\ (f 2 ) be the time position of the previous 
(next) flip (regardless of its type) changing the state of 
any of the bonds bi, ...64 connected to site s. 
Step 4.3 

(T) If ti < t < t2, propose to move the flip to t' = {t% — 
ti) x rndm+ti; 

(2) If t < t 2 < t\ or t 2 < t\ < t, the new time position of 
the flip is t' = mod{(/3 + t 2 — t\) x rndm + ti,/3}. 
Step 4.4 The acceptance ratio is 

R = min{l,cxp{- J AC/(r)dr}} . 
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Update 5 Create/annihilate combinations of one- and 
four-spin flips. 

Step 5.1 Randomly select a site s. 

Step 5.2 With ao kept unchanged, propose to create or 
annihilate a four-spin flip on s with equal probability. 
(T) In creation, draw a random integer k = [(N(s) + 1) x 
rndm] (k = 0, 1 • ■ • N(s)). Propose the time position 

for a new four-spin flip within the interval (t^ , t^ 8 )^ 

using t = (t^j — rjf^) x rndm+T^' 1 . The detailed- 
balance factor associated with this proposal is f c (s) = 

(N(s) + l)(T^ 1 -ri s) )J s . 

(2) In annihilation, exit the update if N(s) = 0; oth- 
erwise draw a random integer k = [N(s) x rndm] + 1 
(k = 1,2- ■■N(s)), and propose to annihilate the kth 
four-spin flip, of which the time position is at t. The 
detailed-balance factor associated with this proposal is 
f a {s) = N(s). Note that in terms of the current config- 
uration parameters r(s) = f c (s)/f a (s) = (t^i — )J S 

in creation and r(s) = f a /f c = l/((r^ ) 1 - T^l\)J s ) in 
annihilation. 

Step 5.3 On bond bi (i=l, 2, 3, 4) connected to site s, 
the time t falls within some time interval (Tj b '\ Tj+l), 
where t is the time position for either the created or the 
annihilated four-spin flip in Step 5.2. For each bond 
bi, propose to create or annihilate a one-spin flip on this 
interval. [As before, «o is kept unchanged, and the prob- 
abilities for creation and annihilation are set equal.] 
(T) For creating a one-spin flip, select new time vari- 
able at random t^O — (t^\ — r j b! ' ) ) x rndm+r,j b '\ 
The corresponding detailed-balance factor is f c (bi) = 

@ For annihilating a one-spin flip, randomly choose be- 
tween labels j and j + 1 and propose to annihilate the 
selected j'th one-spin flip (exit the update if j' = or 
N(bi) + I). The detailed-balance factor is f a (bi) = 2. 
Again, in terms of the current configuration parameters 
r(bi) = fc{bi)/ fa{bi) = (T^l-Tf^)h x /2 in creation and 

r(bi) = 2/((rj, 6 ^ 1 — Ty_\)h x ) in annihilation. 
Step 5.4 The acceptance ratio is 

4 r 

R = xam{l,r(s)Y[r(bi)eiqp{- / A[/(r)dr}} , 

Fig. 6 provides an illustration of configuration changes 
produced in these updates. 

IV. PHASE DIAGRAM 

We restrict our study of Hamiltonian (2) to the pa- 
rameter subspace defined by J s = J p (and use J s as a 
unit of energy). We also set (3 = L in the Monte Carlo 
simulation aimed at the study of ground state properties, 
where L x L = N is the square lattice size. For J s = J p 
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FIG. 6: (Color online) An illustration of Update 5. From 
left to right, the update annihilates a four-spin flip on site s 
and two one-spin flips on bonds b\ and 62 and creates two 
one-spin flips on bonds 63 and 64. 

the location of the self-duality line is given simply by 
h x = h z . Thus, in the following we only need to focus 
on the region h x < h z of the phase diagram. The h x ~0 
case was studied in Ref. 14 where the authors have shown 
that TCM with h z 7^ is equivalent to the 2D transverse- 
field Ising model featuring a continuous phase transition 
with the dynamic critical exponent z = 1. Using high- 
accuracy results for the transvcrsc-ficld Ising model 15 , we 
immediately establish the critical value of the h z -&e\d, 
h { z \h x = 0)^0.328474(3). When h x is finite but weak, 
TCM undergoes the same Ising-type continuous transi- 
tion driven by large /i z -field !> . Though in finite fields the 
predominantly charge-condensed and predominantly vor- 
tex condensed phases can be continuously transformed 
into each other, along the self-duality line (h x = h z ) the 
symmetry between charges and vortices can be broken in 
a first-order transition over a certain range of fields h z 
(=h x ). 

The second-order phase transition line. In a finite h x 
field, our Monte Carlo data show that both (^4 S ) (Fig. 7) 
and {al) (Fig. 8) develop strong non- analytic features at 
the transition point indicating the onset of charge con- 
densation. To compare our data to the series expan- 
sion results, we also compute (A s ) and (a^) using the 
Hellmann-Fcynman theorem for the ground-state energy 
at order 10 in the low-field limit 12 . As can be seen in 
Fig. 7 and Fig. 8, the results of the two approaches match 
very well at low fields, but the series expansion result is 
completely missing the non-analytic behavior in larger 
fields. 

There is no explicit local order parameter character- 
izing the phase transition in terms of original spins 9 ' 10 . 
However, local (in space) order parameter can still be 
found by introducing auxiliary variables and applying 
gauge transformation to the Hamiltonian (2) 9 . As an 
effective magnetic moment associated with operator A s , 
we introduce an order parameter [i s which can be easily 
sampled in the Monte Carlo simulation: 

Ms = ±I{(r< s ) - 0) - (r^ r x W ) + • • • + (-l)^)" 1 
C^-r^ + C-l^WCS-r^,)}, (5) 

(s) (s) (s) 

where r-{ < t 2 < • • • < t n( s } are the imaginary time 
positions of four-spin flips on site s. Since the sign of /z s 
is arbitrary, we define its second and fourth moments as: 

S S 
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The Binder cumulant 1 ' 1 of fx g is denned as 

Q(/i s ) = (^)7</4>- 



(7) 



As shown in Fig. 9, Q(p- S ) behaves linearly in the vicinity 
of the critical point. The finite-size behavior of Q(p. s ) 
near the critical point can be parameterized as 1 ', 



Q (/■*») = Q + Efc=l a k{h z 

+b 1 L^ + b 2 Ly- 2 +dLy s (h z 



h (c)^k L ky t 



(8) 



which includes both the leading scaling terms and cor- 
rections to scaling. The Monte Carlo data are fitted on 
the basis of this formula, according to a least-squared 
criterion. The fitting result shows that the critical point 
for h x = 0.3 is located at h{ c) (h x = 0.3) = 0.333(1). 
Our limited system sizes are insufficient for accurate de- 
termination of the critical exponent yt, which should be 
1.5868(3) 18 since the transition is supposed to be of 3D 
Ising universality. 
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FIG. 7: (Color online) (A s ) versus h z for h x = 0.3. 
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FIG. 8: (Color online) (a§) versus h z for h x =0.3. 

The first-order phase transition line. Our Monte Carlo 
data for (A s ) and (a^) along the self-duality line are re- 
spectively shown in Fig. 10 and Fig. 11, which also dis- 
play the series expansion results both at order 10 in the 
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FIG. 9: (Color online) Binder cumulant Q{n a ) for linear sys- 
tem sizes L = 10, 12, 16, 20 versus h z for h x =0.3. 



low-field limit 12 and at order 5 in the high-field limit 1 ". 
Again, the Monte Carlo data and the series expansion re- 
sult are in good agreement at low fields (approximately 
h z < 0.34). However, the series expansion results do not 
provide an accurate description of the model in the pa- 
rameter region 0.4 < h z < 0.8. As shown in Fig 11, 
the susceptibility grows as the lattice sizes increases at 
intermediate-field region, which is a sign of phase transi- 
tions. 

Along the self-duality line, there is an obvious symme- 
try between charges and vortices. To see if the symmetry 
is broken we sample the probability distribution -P(A) for 
the quantity A = A s — B p . By symmetry -P(A) is sup- 
posed to be an even function of A. The evolution of P{A) 
along the self-duality line is shown in Fig. 12. We observe 
that P (A) at h x = h z =0.3 and 0.6 has a single Gaussian- 
type peak centered at A = 0, while at h x = h z = 0.39 
-P(A) has two peaks at finite values of |A|. The double- 
peak feature of P(A) is a smoking-gun evidence for the 
symmetry breaking between charges and vortices in the 
first-order transition. At h x = h z = 0.42 and 0.43, the 
bimodal structure docs appear for the available lattice 
sizes, but the local minima around A = increases with 
system size and thus should vanish for sufficiently large 
sizes. From finite-size scaling theory, this means that 
despite of the absence of symmetry breaking, these two 
points are close to the lst-order phase transition. To- 
gether with the histogram for h x = h z = 0.41, it is con- 
servative to expect that the end point of the lst-order 
transition line lies in range (0.41,0.42]. 

Symmetry breaking between charges and vortices leads 
to the 2-fold degenerate ground state, which can be clas- 
sified respectively by the predominant condensation of 
charges (A = A s — B p < 0) and predominant conden- 
sation of vortices (A = A s — B p > 0). Since jtx s is used 
as the order parameter for the phase transition to the 
charge-condensed phase we also sample it near the low- 
field end of the first-order transition line but only in con- 
figurations with A = A s — Bp < 0. We use jjL 8 A ^ nota- 
tion to stress that [i s was sampled selectively this way. 
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FIG. 10: (Color online) (A s ) versus h z along the self-duality 
line. 
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FIG. 11: (Color online) (crj) versus h z along the self-duality 
line. The inset shows the susceptibility of al defined as 
XK) = 2L 2 p({(aif) - (atf). 



The behavior of the Binder cumulant Q([i{ A ^) near the 
point where all three transition lines meet is shown in 
Fig. 13. The Monte Carlo data for Q{yi A] ) for differ- 
ent system sizes were also fitted on the basis of Eq. (8), 
showing that the first-order transition line starts around 



,('-) 



0.340(2), in good quantitative agreement 
with the result hi c) = h z c) = 0.3406(4) obtained in 
Ref. 10. The agreement is expected since the two ap- 
proaches give almost the same result for (A s ) and (a^) 
in low field. 

The Binder cumulant of A along the self-duality line is 
shown in Fig. 14. Though Q(A) is highly non-linear near 
the low-field critical point it behaves linearly near the 
high-field end point of the first-order transition line. Near 
this end point, the analysis of intersections of Q(A) for 
different system sizes provides an accurate estimate of the 
critical point (Fig. 15). The extrapolation of intersection 
points for two subsequent sizes is shown in Fig. 16 using 
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FIG. 12: (Color online) Probability distribution of A along 
the self-duality line. -P(A) with bimodal structure are 
reweighted so that the two peaks have equal heights. 
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FIG. 13: (Color online) Q(ji a ) versus h z along the self- 
duality line. 



power-law variable L d for the horizontal axis. Extrap- 
olations for L~ 2 9 , L~ 3 , L~ 31 and L~ 3 2 show that the 
first-order transition line ends at h x = hz^ = 0.418(2). 
This is consistent with our expectation based on his- 
tograms, and differs from the result hie = h z c ^ = 0.48(2) 
of the series expansion 1 " which is questionable in this 
parameter range. It is also possible that part of the dis- 
crepancy is due to systematic error introduced by Dlog 
Padc approximation/extrapolation, used in Ref. 10. Se- 
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FIG. 14: (Color online) Q(A) versus h z along the self-duality 
line. 




0.50 



separate topologically protected ground states from mag- 
netically ordered phases A and B characterized by con- 
densation of charges and vortices, respectively 9,14 . Along 
the blue dashed line in Fig. 17, the symmetry between 
charges and vortices is broken in the first-order transition 
which is accompanied by an abrupt change in the den- 
sity of charges and vortices. However, when large enough 
external field is continuously rotated in the (x, z) plane 
one can continuously convert phases A and B into each 
other without inducing a phase transition. 
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FIG. 15: (Color online) Q(A) versus h z near the end point of 
the lst-order transition line. 



FIG. 17: (Color online) Phase diagram of the toric code model 
in magnetic field. The second-order transition lines are shown 
by full (red) lines and the first-order transition line is repre- 
sented by dashed (blue) line. 



ries expansion to higher order about the high field limit 
and the use of other extrapolation schemes, can provide 
estimates for the end point in agreement with our value 10 . 

The resulting phase diagram of TCM in magnetic field 
is shown in Fig. 17. The continuous phase transition lines 

0.46 | , 1 , 1 , 1 



0.45 - 
0.44 - 
0.43 - 




0.42 

0.000 0.001 0.002 

FIG. 16: (Color online) Extrapolation of intersection points 
for Q(A, L) curves for different system sizes. 



V. SUMMARY 

To simulate the toric code model in a parallel field, 
we developed an algorithm based on the continuous-time 
Monte Carlo approach. It allows one to perform sim- 
ulations of quantum spin models with kinetic/exchange 
terms involving four spins. The phase diagram was con- 
structed using finite-size scaling analysis of Binder cumu- 
lants. The two second-order transition lines enclosing the 
topological phase and related by self-duality were found 
to be in excellent quantitative agreement with the study 
performed in Ref. 10. As for the end point of the first- 
order transition line, our numerical result provides a po- 
tentially better estimate, compared to the result obtained 
in Ref. 10. System size limitations did not allow us to get 
enough resolution to answer the question of how the three 
phase transition lines merge around h z = h x = 0.34 , 
and, in particular, to verify that they all terminate at 
the same point. This longstanding problem requires fur- 
ther investigation. 

Recently, the toric code model has been generalized to 
7Ln degrees of freedom-'". The Monte Carlo algorithm 
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developed in this paper can be used to study the Z3 toric 
code model in a parallel field by adjusting the number 
of flips to be changed in each elementary Monte Carlo 
update. This is an interesting problem to look at. 
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